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Abstract 

The  stabilization  of  high  order  spectral  elements  to  solve  the  transport  equations  for  tracers  in  the 
atmosphere  remains  an  active  topic  of  research  among  atmospheric  modelers.  This  paper  builds 
on  our  previous  work  on  variational  multiscale  stabilization  (VMS)  and  discontinuity  capturing 
(DC)  [Variational  multiscale  stabilization  of  high- order  spectral  elements  for  the  advection- diffusion 
equation  J.  Comput.  Phys.  231  (2012)  7187-7213]  and  shows  the  applicability  of  VMS+DC  to 
realistic  atmospheric  problems  that  involve  physics  coupling  with  phase  change  in  the  simulation 
of  3D  deep  convection.  We  show  that  the  VMS+DC  approach  is  a  robust  technique  that  can  damp 
the  high  order  modes  characterizing  the  spectral  element  solution  of  complex  coupled  transport 
problems.  The  method  has  important  properties  that  techniques  of  more  common  use  often  lack: 
1)  it  is  free  of  a  user-defined  parameter,  2)  it  is  anisotropic  in  that  it  only  acts  along  the  flow 
direction,  3)  it  is  numerically  consistent,  and  4)  it  can  improve  the  monotonicity  of  high-order 
spectral  elements.  The  proposed  method  is  assessed  by  comparing  the  results  against  those  obtained 
with  a  fourth-order  hyper-viscosity  programmed  in  the  same  code.  The  main  conclusion  that  arises 
is  that  tuning  can  be  fully  avoided  without  loss  of  accuracy  if  the  dissipative  scheme  is  properly 
designed.  Finally,  the  cost  of  parallel  communication  is  that  of  a  second  order  operator  which 
means  that  fewer  communications  are  required  by  VMS+DC  than  by  a  hyper-viscosity  method; 
fewer  communications  translate  into  a  faster  and  more  scalable  code,  which  is  of  vital  importance 
as  we  approach  the  exascale  range  of  computing. 
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1.  Introduction 

In  the  field  of  numerical  weather  prediction,  regardless  of  the  numerical  method  that  atmo¬ 
spheric  models  use,  dissipation  of  some  sort  is  added  for  various  reasons.  Stabilization  is  one  of 
them  and  represents  the  main  topic  of  analysis  of  this  paper.  Due  to  its  ease  of  implementation 
and  robustness,  artificial  diffusion  is  the  most  common  mean  of  dissipation  that  is  found  in  cur¬ 
rent  research  and  operational  weather  forecast  models  (cf.  Jablonowski  and  Williamson  (2011) 
and  citations  therein).  Given  an  advection-dominated  problem  described  by  the  general  evolution 
equation 

Qt  =  £q,  (1) 

where  £  is  a  set  of  differential  operators  acting  in  space  on  q,  the  explicit  dissipation  of  q  is  classically 
modeled  as  a  diffusion  operator  (hyper-viscosity,  HV,  from  now  on)  in  the  form  of 

HV  =  (-l)a+1Va-(u2aVaq),  (2) 

where  a  is  a  positive  integer  and  a  is  the  diffusivity  coefficient.  HV  is  sometimes  found  also 
in  other  fields  of  computational  fluid  dynamics;  one  example  is  found  in  Baruzzi  et  al.  (1992)  to 
preserve  higher  order  accuracy  in  the  simulation  of  high  speed  flows.  In  the  stabilization  of  finite 
elements,  the  orthogonal  subscale  (OSS)  approach  for  sub-gridscale  modeling  (see,  e.g.,  Codina 
(2002))  could,  in  principle,  also  be  seen  as  a  dynamic  version  of  such  a  scheme.  However,  even  if 
HV  in  the  form  of  Eq.  (2)  is  scale-selective  and  only  damps  higher  frequencies,  it  is  not  physical. 
With  respect  to  consistency,  the  diffusion  coefficient  should  decrease  as  the  grid  is  refined;  this 
condition  is  respected  by  HV  if  the  value  of  its  coefficients  is  reduced  as  the  resolution  increases, 
approaching  zero  as  the  grid  is  sufficiently  refined.  Nonetheless,  because  HV  does  not  depend  on 
the  residual  of  the  equations,  its  strength  and  anisotropy  cannot  adapt  to  the  flow  direction  and 
strength  automatically.  Additionally,  for  a  stabilizing  scheme  to  preserve  the  shape  of  the  tracer, 
dissipation  normal  to  the  flow  should  be  avoided  whereas  the  dissipation  in  the  direction  of  the  flow 
should  be  sufficient  to  control  instabilities  or  unboundedness  (Hughes  and  Brooks,  1979;  Hughes 
et  ah,  1986;  Codina,  1993).  This  condition  will  not  be  respected  by  diffusion  schemes  such  as  HV, 
as  exemplified  in  Fig.  1.  To  preserve  the  correct  physical  dimensions  of  the  hyper- viscous  term, 
the  value  of  U2a  must  scale  correctly  with  respect  to  a  and  the  grid  spacing.  Its  selection  not  only 
is  non-trivial,  but  has  a  great  impact  on  the  solution  of  the  problem.  Quoting  Jablonowski  and 
Williamson  (2011)  ’’The  choice  of  the  V2,  V4  or  even  higher-order  diffusion  coefficient  is  most  often 
motivated  by  empirical  arguments  and  chosen  in  a  somewhat  arbitrary  manner.  It  is  sometimes 
even  considered  a  model  tuning  parameter  [...]”. 

It  is  well  known  that  the  straightforward  numerical  approximation  of  Eq.  (1)  that  models 
advection-dominated  problems  is  characterized  by  Gibbs  oscillations  that  may  render  the  solution 
unacceptable.  These  oscillations  are  treated  by  either  anti-aliasing  filters  (  e.g.,  Vandeven  (1991); 
Boyd  (1998))  alone  or,  possibly,  by  a  combination  of  the  latter  with  some  type  of  artificial  diffusion 
in  the  form  of  (2)  to  improve  the  non-monotonic  filtered  solution.  Research  in  finding  the  optimal 
solution  to  this  problem  is  still  ongoing  (see,  e.g.,  Malm  et  al.  (2013)). 
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(Min.Max)  =  (0.0.  0.9482) 


(Mln.Max)  =  (-0.0739.  1.1052) 


(Mill.  Max)  =  (-4.125e-06.  1) 


Figure  1:  Advection  of  a  2D  square  wave  in  a  doubly-periodic  channel  with  velocity  directed  along  the  x-axis,  which 
corresponds  to  the  bottom-right  edge  of  the  squared  surfaces  in  the  three  plots.  Stabilization  achieved  by  the  operators 
V2  (left)  and  V4  (center),  with  constant  coefficients  V2  =  0.001  m2  s_1  and  1/4  =  0.000001  m4  s_1,  respectively,  and 
using  VMS  (right).  Reprinted  from  Marras  et  al.  (2012),  Fig.  21,  with  permission  from  Elsevier. 


In  the  quest  for  a  solution  to  the  drawbacks  of  HV  and  for  the  improvement  of  the  monotonicity¬ 
preserving  properties  of  stabilized  high-order  methods,  the  main  novelty  of  our  current  paper  lies 
in  the  application  of  the  variational  multiscale  stabilization  (VMS)  by  Hughes  (1995)  with  the 
crosswind  discontinuity  capturing  by  Codina  (1993)  (DC  from  now  on)  to  stabilize  the  spectral 
element  solution  of  three  coupled  advection  equations  in  the  simulation  of  fully  3D  deep  convection. 
To  our  knowledge,  the  only  application  of  VMS  to  high  order  spectral  elements  for  coupled  problems 
appears  in  Gjesdal  et  al.  (2009)  for  the  turbulence  modeling  of  incompressible  flows.  In  the  case  of 
Marras  et  al.  (2012),  VMS+DC  was  only  applied  to  solve  one  scalar  advection  problem  and  verified 
against  standard  benchmarks  in  one  and  two  dimensions.  The  characteristic  velocities  of  deep 
convection  may  be  orders  of  magnitude  larger  than  those  in  the  purely  academic  tests  presented 
in  our  previous  work,  thereby  subjecting  the  scheme  to  a  more  demanding  test  when  it  comes  to 
stabilization  and  treatment  of  over-  and  undershoots.  In  the  simulation  of  deep  convection,  this 
scheme  was  successfully  applied  by  Marras  et  al.  (2013)  for  2D  problems  using  linear  finite  elements. 
The  main  goal  of  this  work  is  then  to  prove  the  applicability  of  VMS+DC  to  3D  realistic  atmospheric 
simulations  and  show  that  VMS+DC  is  1)  parameter-free,  2)  anisotropic,  3)  numerically  consistent, 
4)  inexpensive  on  parallel  computers,  and  5)  can  improve  the  monotonicity-preserving  properties 
of  high  order  spectral  elements.  We  implement  this  stabilization  method  within  the  Nonhydrostatic 
Unified  Model  of  the  Atmosphere  (NUMA)  (Kelly  and  Giraldo,  2012;  Giraldo  et  al.,  2013). 

The  equation  sets  and  the  details  of  the  residual-based  stabilization  and  crosswind  dissipation 
are  reported  in  Sec.  2.  The  fundamental  results  for  a  fully  3D  squall  line  simulation  are  shown  in 
Sec.  3.  A  summary  of  our  conclusions  is  discussed  in  Sec.  4. 

2.  Equations,  their  discretization,  and  stabilization  techniques 

Let  SI  6  i3  be  a  fixed  three  dimensional  domain  with  boundary  dUl  and  Cartesian  coordinates 
x  =  (x,  y,  z).  Let  us  identify  the  dry  air  density,  the  velocity  vector,  and  the  potential  temperature 
with  the  symbols  p,  u,  and  6.  Let  us  also  define  the  mixing  ratios  of  water  vapor,  cloud  water, 
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and  rain  as  qv  =  pv/p,qc  =  Pc/ P  and  qr  =  pr/p,  where  pv,c.r  are  the  densities  of  water  vapor, 
cloud,  and  rain.  Furthermore,  let  p'(t,x)  =  p(t,x)  —  po{z),  =  9(t,x)  —  9o(z),  and  p'(t,x )  = 

p(t,x )  —  po(z)  be  the  perturbations  of  density,  potential  temperature,  and  pressure  with  respect 
to  a  hydrostatically  balanced  background  state  indicated  by  the  subscript  0  (see,  e.g.,  Marchuk 
(1974)).  Then,  the  strong  form  of  the  time-dependent  Euler  equations  that  model  a  stratified, 
moist  atmosphere  can  be  written  as1: 

p't  +  u  •  Vp  +  pV  •  u  =  0, 


ut  +  u-  Vu  +  f V  •  (Ip') 


0k(l  +  eqv  -qc-  qr), 


9't  +  u-Vd 


Se(p,  9,  qv,  qc,  qr), 


(3) 


qit  +  u-\7qi  =  Sqi  (p,  9,  qv,  qc,  qr),  for  i  =  v,c,r, 

where  I  is  the  rank-3  identity  matrix,  gh  is  the  acceleration  of  gravity  directed  as  k  =  [0  0  l]7”, 
and  e  =  R/Rv  is  the  ratio  of  the  gas  constants  of  dry  air,  R,  and  of  water  vapor,  Rv.  9 ,  p,  and 
p  are  related  through  the  equation  of  state  for  a  perfect  gas.  Because  moist  air  contributes  to 
the  buoyancy  of  the  flow,  the  right  hand  side  of  the  momentum  equation  is  corrected  by  the  total 
buoyancy  B  =  gk(l  +  eqv  —  qc  —  qr)-  Due  to  the  nricrophysical  processes  that  involve  phase  change 
in  the  water  content,  the  source/sink  terms  identified  by  Sg!gi  are  modeled  and  computed  via  the 
Kessler  microphysics  for  warm  clouds  (Kessler,  1969).  The  system  of  Eqs.  (3)  must  be  solved  in  S7 
Vi  G  (0 ,tf)  with  properly  assigned  initial  and  boundary  conditions.  Because  in  this  study  we  are 
mainly  interested  in  the  solution  of  the  advection  equations  of  water  tracers,  we  concentrate  on  the 
last  three  equations  of  the  system.  To  simplify  the  description  of  their  numerical  treatment  later, 
we  re-express  them  in  compact  form  as: 


qt  +  £q  =  S(q), 


(4) 


where 


Qv 

Qv 

'V 

qc 

,  £q  =  u  •  V 

Qc 

t  s(q)  = 

Sq  c 

<lr_ 

_qr_ 

[sq J 

(5) 


Because  the  coupling  of  (4)  with  the  Euler  equations  is  done  by  the  saturation  adjustment  of 
Soong  and  Ogura  (1973),  the  solution  is  first  obtained  for  the  homogeneous  counterpart  of  Eq. 
(4)  before  the  sources  are  computed  and  the  unknowns  updated  through  Kessler.  Because  we 


1For  the  sake  of  precision,  the  Euler  equations  are  only  the  first  three  equations  of  the  system,  and  are  coupled  to 
the  equations  for  the  quantities  g;  through  velocity  and  the  source  terms. 
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already  describe  the  spectral  element  approximation  of  the  scalar  counterpart  of  Eq.  (4)  in  our 
previous  work  (Marras  et  al.  (2012)),  we  can  omit  the  details  on  spectral  elements  and  proceed  to 
the  description  of  their  stabilization. 

2.1.  Stabilization  of  spectral  elements 

Let  us  define  the  space  of  trial  solutions  as  Q  =  {q  €  H1(fl)  :  q  =  (/on <914}  and  the  space  of 
test  functions  =  {if  €  if1  (14)  :  if  =  Oon<914},  where  bold  indicates  a  vector  quantity  and  the 
discrete  counterpart  of  Q  and  4'  will  be  indicated  by  the  superscript  h.  Given  the  basis  function 
ifh,  the  stabilized  spectral  element  solution  of  (4)  consists  in  finding  the  solution  q/(  £  Qh  such 
that 


r 


qf+Cqh  <mh  +  b(iph,qh) 


iphS(q  )dnh 


holds  V  iph  e  yih.  In  Eq.  (6), 


(6) 


h,qh) 


Ne 

— E 

e=l  ' 


C*iph  TU(qh)dnh 


is  the  stabilizing  term  defined  as  a  function  of:  the  adjoint  operator  C*  acting  onto  the  basis 
function,  the  equations  residual  ^(q^-),  and  the  parameter 


T  = 


2|u| 


-1 


+  epsilon 


where  h  is  a  characteristic  length  of  the  computational  grid.  Because  of  the  somewhat  uncertain 
definition  of  r  (Knobloch  (2008)),  the  choice  of  h  may  vary  as  a  function  of  the  physical  characteris¬ 
tics  of  the  problem  and  of  the  structure  of  the  grid.  For  example,  Owen  et  al.  (2013)  achieved  better 
solutions  by  taking  h  as  the  length  of  the  shortest  element  edge  in  the  flow  boundary  layer.  In  this 
study  we  use  h  =  \J Ax2  +  Ay2  +  A z2,  where  A-  is  the  mean  distance  between  two  consecutive 
nodes  within  a  high  order  element  with  unequally  spaced  nodes.  Because  velocity  may  be  zero,  the 
machine  epsilon  is  added  to  avoid  divisions  by  zero  on  the  computer.  Should  it  be  the  case  that 
Eq.  (4)  included  some  physical  diffusion  v  (e.g.  eddy  viscosity,  thermal  diffusion),  r  would  simply 
include  it  in  the  form 


v  2 1  u|  \ 

T  =  1  W  +  IT  +  epsil°”) 


This  is  not  the  only  definition  of  r  that  can  be  found  (see,  e.g.,  Knobloch  (2008)).  For  up  to 
cubic  elements  with  equally  spaced  nodes,  a  derivation  of  the  components  of  r  was  obtained  by 
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Table  1:  Recommended  values  of  the  constant  C  as  a  function  of  the  elements’  order, 
linear  quadratic  cubic  quartic 

~C  07  035  023  0.175 


Houzeaux  et  al.  (2009).  Based  on  Houzeaux  et  al.,  r  for  elements  with  unequally-spaced  nodes  was 
obtained  by  Marras  et  al.  (2012)  although,  for  the  sake  of  simplicity  and  without  a  great  loss  of 
accuracy,  it  is  not  used  in  this  study.  Because  b{t^h ,  qft  )  is  a  function  of  the  equation  residual  and  of 
a  characteristic  grid  size,  the  effect  of  stabilization  goes  to  zero  as  the  exact  solution  is  approached. 
This  makes  the  method  numerically  consistent.  One  additional  property  of  this  stabilizing  scheme 
is  its  anisotropy.  The  stabilizing  effect  only  acts  along  the  direction  of  the  flow,  which  is  the  only 
direction  where  advection  dominates  diffusion  (if  any). 


2.2.  Crosswind  discontinuity  capturing 

Stabilization  in  the  form  of  (6)  is  neither  monotonic  nor  monotonicity  preserving;  hence,  over- 
and  undershoots  are  still  likely  to  affect  the  solution.  This  is  especially  true  for  large  gradients  in  the 
transported  quantities.  To  overcome  this  issue  for  low  order  finite  elements,  Codina  (1993)  designed 
an  additional  controlled  crosswind  discontinuity  capturing  that  acts  in  the  proximity  of  internal 
and  boundary  layers  and  that  is  orthogonal  to  the  wind  direction.  The  additional  dissipation  to  be 
added  to  Eq.  (6)  is  built  as 


Ne 


DC  =  Y. 


e=l 


[  1 

/  -max 

lnh  2 


0  ,C- 


2v 


u||l  hj  |V<?' 


U  (X)'  u 


U 


Vqh  dnh 


(7) 


where  U||  is  the  velocity  component  in  the  streamline  direction  and  <g>  indicates  a  tensor  product. 
The  constant  C  is  not  user-defined.  Its  values  were  suggested  by  Codina  (1993)  for  linear  and 
quadratic  elements.  Based  on  the  values  provided  by  Codina,  we  extend  them  to  higher  order 
polynomials  by  a  simple  extrapolation  from  the  linear  and  quadratic  values  and  report  them  in 
Table  1.  Note  that  Codina  assumes  finite  elements  with  equi-spaced  nodes.  Although  the  values 
that  he  suggests  seem  to  conform  sufficiently  well  with  the  Legendre-Gauss-Lobatto  points  at  order 
4,  a  proper  study  should  be  made  to  define  C  for  higher  order  spectral  elements.  To  our  knowledge, 
this  remains  an  open  topic  but  falls  beyond  the  scope  of  this  paper.  From  now  on,  the  symbol 
VMS+DC  will  be  used  to  indicate  stabilization  with  crosswind  discontinuity  capturing. 


2.3.  Anisotropic  HV 

Because  of  the  large  aspect  ratio  of  the  grid  elements  that  we  use  in  the  three-dimensional 
cloud  simulations  presented  below,  where  hx,y  3>  hz,  we  implemented  an  anisotropic  version  of  the 
hyper- viscosity  operator  defined  in  Eq.  (2)  and  write  it  as 

HVa  =  (-1)“+1VQ  •  (V2avaq),  (8) 
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where 


0 


0 

0 


v  = 


+: 

0  Vy 

Q  0  vz 


is  the  matrix  of  the  diffusion  coefficients  along  the  x,  y  and  z  directions.  This  scheme  may  be 
seen  as  a  simplified  version  of  the  tensor  hyper-viscosity  that  was  recently  analyzed  in  Guba  et  al. 
(2014).  The  implementation,  albeit  simple,  requires  some  attention  if  the  higher  order  operator  is 
built  recursively.  After  building  the  second  order  Laplace  operator  using  Eq.  (8),  v  must  be  set  to 
the  identity  matrix  for  every  subsequent  application  of  the  higher  order  derivation. 

Remark  1.  The  anisotropy  of  HVa  must  be  understood  in  a  different  way  with  respect  to  the 
anisotropy  of  VMS+DC.  HVa  owes  its  anisotropy  to  the  fact  that  the  magnitude  of  the  dissipation 
varies  based  on  the  shape  of  the  grid  but  it  still  acts  in  all  directions  regardless  of  the  flow  features. 
On  the  other  hand,  VMS+DC  is  anisotropic  in  that  the  dissipation  acts  exclusively  along  the  flow 
direction  and  not  normal  to  it.  This  is  an  important  feature  that  stabilizing  schemes,  unlike  physical 
stresses  or  diffusion,  should  preserve  to  least  affect  the  flow  physics. 


3.  Numerical  experiments 

To  verify  that  the  NUMA  implementation  of  the  method  is  correct  and  able  to  reproduce  the 
results  of  our  previous  work,  we  first  run  the  Atmo-3D  benchmark  of  Marras  et  al.  (2012).  The 
simulation  of  a  fully  three-dimensional  squall  line  follows  Atmo-SD.  To  understand  how  the  method 
performs  on  moist  dynamics,  the  negative  values  of  the  three  moisture  variables  are  not  filtered  in 
the  Kessler  microphysical  scheme  as  is  typically  done  to  avoid  negative  moisture. 


3.1.  Atmo-3D:  cylindrical  tracer  in  a  buoyant  atmosphere 

A  passive  tracer  is  transported  in  a  neutrally  stratified  atmosphere  contained  in  12  =  [1  x  oo  x  1] 
km3.  The  simulation  final  time  is  tf  =  600  s.  The  neutral  background  atmosphere  at  uniform 
potential  temperature  6q  =  300  K  is  perturbed  by  a  cylindrical  thermal  of  radius  rc  =  250  m, 
centered  in  (xc,  zc)  =  (500,  350)  m,  aligned  with  the  y-axis,  and  defined  by 


6'  = 


0.5 

0 


1  +  cos 


if  r  <  rc 
elsewhere, 


(9) 


where  r  =  \/{x  —  xc )2  +  (z  —  zc)2.  The  east,  west,  top,  and  bottom  boundaries  are  modeled  as 
free-slip  solid  walls  whereas  periodic  boundary  conditions  are  imposed  on  the  north  and  south  walls 
(y-direction) . 

The  results  plotted  in  Fig.  2  are  quantitatively  and  qualitatively  comparable  to  the  simulations 
presented  in  our  previous  work  where  the  element  grids  are  defined  by  10x1x10  and  20x1x20 
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Figure  2:  Atmo-3D:  passive  tracer  in  a  rising  thermal  field.  Vertical  slice  of  qtracer  at  y  =  0  m.  Left:  initial  condition. 
Center:  coarse  simulation  with  10  x  1  x  10  elements  of  order  4.  Right:  fine  simulation  with  20  x  1  x  20  elements  of 
order  4.  The  results  are  shown  at  time  tf  =600  seconds. 


elements  of  order  4.  Nonetheless,  a  careful  eye  may  spot  minor  differences.  The  differences  are 
due  to  the  use  of  different  time-integrators  between  the  two  codes,  which  implies  different  time 
steps.  The  velocity  held  that  drives  the  passive  tracer  is  obtained  by  the  coupling  of  the  advection 
equations  with  the  Euler  equations  that  model  the  motion  of  the  thermal  bubble.  To  preserve  a 
numerically  stable  SEM  solution  of  the  Euler  equations,  the  Boyd-Vandeven  filter  (Vandeven,  1991; 
Boyd,  1998)  was  applied.  Because  this  filter  is  not  idempotent,  its  application  at  different  time 
steps  will  give  different  solutions  of  thermal  perturbation  and,  consequently,  affect  the  solution  of 
the  water  tracers.  However,  the  results  between  the  two  model  simulations  are  sufficiently  similar 
which  gives  us  confidence  that  the  implementation  of  VMS+DC  in  NUMA  is  coded  properly.  We 
would  like  to  stress  that  no  filter  was  applied  to  the  spectral  element  solution  of  the  advection 
equations  for  the  three  tracers,  that  were  solely  stabilized  by  VMS  and  DC.  We  now  proceed  and 
verify  the  methodology  in  a  fully  coupled  and  convection-permitting  system  where  the  tracers  are 
subjected  to  phase  change  and  the  equations  are  thermodynamically  coupled  to  the  Euler  equations 
of  compressible  flows. 

3.2.  Squall  line 

To  further  evaluate  VMS+DC  we  apply  it  to  the  simulation  of  (I)  the  squall  line  of  Weisman 
et  al.  (1988)  (WKR88,  from  now  on)  and  of  (II)  the  squall  line  over  a  flat  terrain  of  Frame  and 
Markowski  (2006)  (FM06,  from  now  on).  The  two  simulations  are  initialized  with  the  same  sounding 
given  by  Weisman  and  Klemp  (1982)  but  differ  in  the  domain  size  and  in  the  intensity  of  the  thermal 
perturbation.  A  westerly  sheared  flow  has  an  increasing  velocity  ranging  from  0  ms-1  at  the  surface 
to  17.5  ms-1  at  z  —  2.5  km  and  upward.  We  consider  it  useful  to  execute  both  simulations  to 
test  the  sensitivity  of  the  parameter-free  scheme  in  different,  yet  similar,  situations.  Convection  is 
triggered  by  a  cylindrical  thermal  aligned  with  the  y-axis  and  defined  as 
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o'  = 


9C  cos  (^y)2  +  0.5  |sin  (^)  +  cos(y2a3)| 

0 


if  r  <  1, 
elsewhere, 


(10) 


where  r  =  ^ (x  —  xc)2/r%  +  (z  —  zc)2/r%,  9C  =  2  K  in  WKR88,  and  9C  =  4  K  in  FM06.  In  both 
cases,  the  thermal  has  radii  ( rx ,  rz)  =  (10  km,  1.4  km)  and  is  centered  in  the  middle  of  the  domain  at 
(xc,  yc,  zc )•  The  last  two  terms  on  the  right  hand  side  of  Eq.  (10)  define  a  non-periodic  perturbation 
to  the  otherwise  axisymmetric  thermal.  We  use  this  expression  instead  of  a  random  perturbation  of 
9  for  the  sake  of  reproducibility  of  the  same  initial  state.  At  the  four  lateral  boundaries  periodicity 
is  used,  whereas  the  ground  and  top  boundaries  are  modeled  as  inviscid  solid  surfaces.  To  absorb 
the  vertically  propagating  gravity  waves  we  add  a  5  km  absorbing  sponge  at  z  =  zs  =  12  km. 
The  sponge  is  a  Rayleigh-like  layer  defined  as  follows.  In  the  area  occupied  by  the  sponge,  the 
solution  variables  q  are  corrected  as  q  =  T>(q  —  q^),  where  b  indicates  the  value  of  q  on  the  physical 
boundary  and  T>  is  the  damping  coefficient  given  by 


f  sin4 

7 r  (z-zs) 

2  ( Ztop  Zs) 

l  o 

if  z  >  zs 
\1  z  <  zs. 


3.2.1.  Squall  line  I:  WKR88 

For  the  WKR88  test,  we  divide  the  domain  D=[180xl20xl7.5]  km3  into  36x16x20  elements 
of  order  4.  The  simulation  is  executed  until  tf  =14400  seconds.  Without  coefficients  to  be  selected, 
only  one  simulation  is  executed  using  VMS+DC.  On  the  other  hand,  four  simulations  are  run  using 
hyper-diffusion:  a  first  set  of  three  using  HV  with  uniform  coefficients  zq  =  { 104 , 106, 108}  m4s~4 
and  one  using  HVa  with  coefficients  vx  =  vy  =  1  X  10'  m4s_1,z q  =  3  x  104m4s_1.  As  we  will 
notice  below,  these  values  of  zq,  vy  and  uz  would  make  the  simulation  unstable  if  used  isotropically. 
The  time  series  of  maximum  and  minimum  vertical  velocities  are  plotted  in  Fig.  3.  For  all  the 
diffusion  configurations  the  first  updraft  is  triggered  within  the  first  30  minutes.  In  the  case  of 
VMS+DC,  peak  values  between  35  and  45  ms-1  are  reasonably  comparable  to  the  observed  30 
ms-1  of  Weisman  et  al.  (1988)  for  the  same  17.5  ms^1  shear.  The  absolute  value  and  trend  are 
also  in  agreement  with  Weisman  and  Klemp  (1982)  where  a  sequence  of  storm  growths  and  declines 
are  observed  due  to  the  continuous  cycle  of  cloud  formation,  rainout,  and  dissipation.  A  steadily 
decreasing  intensity  of  the  updrafts  is  predicted  by  all  the  HV  solutions  as  well,  although  the 
maximum  strength  of  the  updrafts  larger  than  60  ms-1  is  possibly  indicating  that  none  of  the  HV 
coefficients  selected  so  far  is  sufficiently  dissipative  under  these  flow  conditions.  Albeit  educated,  the 
choice  of  zq  is  practically  arbitrary.  Of  these  three  values,  the  only  zq  that  preserves  stability  for  the 
full  length  of  the  14400  second  simulation  is  zq  =  1  x  106  m4 s_1,  even  if  the  large  85  ms-1  peak  at 
approximately  2.5  hours  and  another  large  peak  at  almost  4  hours  are  the  symptoms  of  an  incipient 
breaking  that  was  eventually  controlled  by  the  end  of  the  simulation.  Given  the  current  grid,  the 
other  values  of  zq  are  either  too  small  or  too  large  and  hence  cause  the  simulation  to  lose  stability. 
In  simple  words,  when  zq  is  too  small,  the  stability  loss  is  due  to  the  insufficient  dissipation  at  the 
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Figure  3:  WKR88:  time  evolution  of  max(w)  (top)  and  min(«;)  (bottom).  Because  the  HV  simulation  with  = 
1  x  10s  m4s-1  lost  stability  within  the  very  first  seconds,  its  time  series  are  not  visible  in  the  plots.  This  is  also 
reflected  in  the  curves  of  qc  and  qr  in  Fig.  4. 
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given  resolution.  The  value  z^4  =  1  x  104m4s-1  preserves  a  stable  simulation  for  approximately  2 
hours  but  breaks  right  after,  as  visible  from  the  trimmed  black  curve  in  Fig.  3.  The  problem  can  be 
solved  by  either  increasing  the  resolution  or  increasing  the  dissipation  coefficient.  In  either  case,  the 
time  step  must  be  decreased  to  respect  the  stringent  limit  that  diffusion  and  resolution  impose  on 
the  time  step.  Although  we  are  using  an  implicit-explicit  (IMEX)  time  integration  scheme  to  solve 
the  system  of  equations  (3),  diffusion  is  currently  treated  explicitly.  Diffusion  may  be  included  in 
the  IMEX  time-integrators  but  it  would  prevent  us  from  deriving  a  Schur  form,  which  is  the  reason 
for  the  success  of  IMEX  methods.  Because  it  falls  beyond  the  scope  of  this  paper,  we  will  not  enter 
the  details  of  time  integration  and  its  analysis;  the  interested  reader  should  refer  to  Giraldo  et  al. 
(2013)  for  more  on  this  topic. 

We  are  aware  that  we  cannot  increase  z^4  arbitrarily  in  this  simulation  because  of  the  highly 
anisotropic  grid  made  of  elements  that  are  much  smaller  along  the  z-direction  than  they  are  along 
x  and  y.  Based  on  scaling  arguments,  this  simply  translates  to  choosing  a  smaller  viscosity  along 
z  but  still  satisfying  the  larger  dissipation  needs  in  the  two  other  directions.  We  then  execute  the 
simulation  using  the  HVa  (Eq.  (8))  with  coefficients  ux  =  vy  =  1  X  10'm4s_1,t'z  =  3  X  104m4s_1 
to  reach  a  stable  solution  that  admits  maximum  updrafts  within  the  acceptable  limits  presented  in 
previous  studies.  Improved  results  obtained  with  HVa  with  respect  to  HV  are  noticeable  by  looking 
at  the  the  red,  dashed-dotted  curve  in  Fig.  3. 

With  respect  to  the  treatment  of  the  undershoots,  Fig.  4  shows  that  VMS+DC  greatly  reduces 
the  undershoots  that  affect  cloud  water  and  rain  and  that  are  not  controlled  by  any  of  the  4ih-order 
hyper- viscosities.  To  verify  that  the  curtailing  of  the  undershoots  is  not  the  mere  result  of  excessive 
dissipation  by  VMS+DC,  we  plot  the  evolution  of  the  total  accumulated  rainfall  at  3600  s,  4800  s, 
and  7200  s  in  Fig.  5.  Because  Weisman  et  al.  (1988)  do  not  report  the  values  of  either  accumulated 
or  instantaneous  precipitation,  we  refer  to  the  values  of  total  accumulated  rainfall  reported  by 
Frame  and  Markowski  (2006).  In  spite  of  the  different  domain  and  resolution,  the  values  of  Frame 
and  Markowski  are  still  a  valuable  source  of  expected  precipitation  given  a  certain  initial  sounding. 
In  Frame  and  Markowski  (2006),  after  18000  s  the  precipitation  ranges  between  5  and  50  mm.  In 
this  study,  at  7200  s  the  accumulated  rainfall  ranges  between  0  and  35  mm  (Fig.  5)  and  reaches  a 
maximum  of  60  nun  after  14400  s  (not  shown).  We  should  note  that  by  14400  s  our  squall  line  has 
crossed  the  western  boundary  and  has  re-entered  the  domain  through  the  eastern  boundary  due 
to  periodicity.  For  this  reason,  a  direct  comparison  may  be  no  longer  possible  since,  at  the  second 
pass  of  the  storm,  more  rain  has  precipitated  in  a  region  where  it  would  have  not  precipitated 
otherwise.  The  contours  of  the  accumulated  rain  from  the  HVa  simulation  are  plotted  in  the  same 
figure.  As  partially  expected  from  the  observations  of  the  time  evolution  of  velocity,  qc,  and  qr,  the 
extremely  small  dissipative  effect  of  HVa  causes  a  lot  more  rain  to  form  and  precipitate.  Again, 
this  is  classically  solved  by  better  tuning  the  coefficients. 

In  Fig.  6,  we  show  a  vertical  cross  section  of  the  line-averaged  cloud  content  resolved  by 
VMS+DC  and  HV.  In  agreement  with  the  time  series,  as  time  evolves  from  1800  to  7200  seconds 
the  VMS+DC  solution  preserves  stability  and  smoothness  of  the  cloud  water  content.  The  classical 
anvil  shape  of  the  cloud  is  preserved  by  both  simulations,  although  the  HV  solution  presents  an 
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Figure  4:  WKR88:  time  evolution  of  min((/c)  (top)  and  min(gr)  (bottom).  Because  the  HV  simulation  with  v 4  = 
1  x  108  m4  s'1  lost  stability  within  the  very  first  seconds,  its  time  series  is  not  visible  in  the  plots. 
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Figure  5:  WKR88:  total  accumulated  rainfall  (mm)  at  t/  =  3600s  (top  row),  tf  =  4800s  (middle  row),  and  t/  =  7200s 
(bottom  row).  Left  column:  VMS+DC.  Right  column:  HVa  with  vx  —  vy  =  1  x  10'  m4s~1,i'z  =  3  x  104m4s-1. 
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overall  increasing  oscillatory  behavior  that  eventually  breaks  the  anvil  and  the  overall  shape  of  the 
cloud.  This  is  simply  due  to  the  suboptimal  choice  of  the  diffusion  coefficient.  The  same  occurs 
with  HVa  as  is  visible  in  the  time  series.  The  large  extrema  in  the  time  series  indicate  that  this 
particular  HY  simulation  is  close  to  diverging  and  that  a  slightly  larger  coefficient  will  improve  the 
solution. 

Remark  2.  Throughout  this  article  we  often  give  hints  on  the  possible  solutions  that  could  im¬ 
prove  either  HY  or  HVa,  but  we  limit  ourselves  to  find  the  first  possible  coefficient  that  avoid  the 
simulation  to  diverge.  We  never  go  further  than  this  because  the  purpose  of  this  article  is  not 
that  of  assessing  HV,  as  this  has  been  done  thoroughly  by  other  authors  before.  However,  we 
still  consider  it  necessary  to  show  the  results  obtained  with  HV  and  HVa  to  have  a  direct  form  of 
comparison  of  what  can  be  achieved  with  VMS+DC  when  HV  is  coded  in  the  same  exact  code  and 
the  same  simulations  are  run.  No  conclusive  remarks  on  HV  are  meant  to  arise  from  this  analysis. 

In  Fig.  7  we  are  looking  at  the  time  and  space  evolution  of  the  3D  iso-surfaces  of  qc.  The 
geometry,  the  position  of  the  lifting  condensation  level,  and  the  vertical  extension  of  qc  obtained  with 
VMS+DC  are  in  good  agreement  with  Weisman  et  al.  (1988)  and  Weisman  and  Rotunno  (2004). 
Since  an  analytic  solution  to  this  problem  does  not  exist,  we  cannot  draw  definitive  conclusions; 
however,  given  the  specific  initial  sounding,  certain  properties  of  the  space  distribution  of  the  storm 
should  be  consistent  with  previous  studies,  and  we  have  been  able  to  achieve  this  with  VMS+DC.  A 
first  conclusion  can  already  be  drawn  at  this  stage:  without  a  constant  to  be  tuned,  VMS+DC  can 
still  produce  suitable  solutions.  If  we  count  the  number  of  simulations  that  was  necessary  to  run 
before  obtaining  a  meaningful  result  from  a  parameter-dependent  method,  the  enormous  advantage 
of  a  parameter- free  scheme  is  evident. 

Given  the  same  initial  field  but  a  different  domain  and  different  resolution,  we  then  proceed 
with  the  FM06  test  for  further  testing  that  may  confirm  this  statement. 

3.2.2.  Squall  line  II:  FM06 

For  the  FM06  test,  we  divide  the  domain  D=[300x60x  17.5]  km3  into  80x12x9  elements  of 
order  4.  Our  domain  extends  100  km  less  than  the  FM06  along  the  x-axis.  This  makes  the 
simulation  relatively  smaller  but  without  compromising  the  full  development  of  the  physical  features 
of  the  squall  line.  Periodic  boundary  conditions  are  used  at  the  four  lateral  boundaries  (unlike  in 
Frame  and  Markowski  (2006),  where  periodicity  is  imposed  along  y  only).  The  simulation  is  run 
until  tf  =18000  seconds.  We  execute  the  same  simulation  using  VMS+DC,  isotropic  HV  with 
tq  =  { 104 ,  106, 108}  m4  s'1,  and  HVa  with  vx  =  vy  =  109,  vz  =  3  x  105  m4  s'1. 

The  initial  convection  forms  within  the  first  900  s  for  all  cases;  a  line  of  single  convective  cells 
has  formed  by  1800  s  and  continues  to  develop  while  merging  into  one  single  line  as  the  simulation 
continues  (Fig.  8).  The  same  features  of  the  evolving  storm  are  observed  for  VMS+DC  and  for 
every  stable  HV,  including  HVa.  This  is  stressed  by  the  time  series  of  the  maximum  and  minimum 
vertical  velocities  plotted  in  Fig.  9.  Except  for  1/4  =  1  x  104m4s_1  and  q  =  1  x  106m4s_1,  the 
maximum  vertical  velocities  obtained  with  VMS+DC,  HV  with  1/4  =  1  x  108  m4  s'1,  and  HVa  show 


14 


Z  (m)  Z  (m)  Z  (m) 


Qc,ave  (g/kg).  7200s 


qc,ave  {g/kg).  3600s 


2 


1 


Qc,ave  (g/kg).  7200s 


Figure  6:  WKR88:  line  averaged  vertical  cross  sections  of  qc  at  1800,  3600,  and  7200  seconds.  Stabilization  achieved 
by  VMS+DC  (left  column)  and  by  HV  with  i/4  =  1  x  106m4s_:L  (right  column). 
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Figure  7:  WKR88:  perspective  view  of  the  squall  line  (grey  scale)  at  1800,  3600,  and  7200  seconds.  Stabilization 
achieved  by  VMS+DC  (left  column)  and  by  HV  with  v 4  =  lx  10®  m4s-1  (right  column).  The  velocity  vectors  at  the 
surface  are  plotted  every  7  grid  points.  The  instantaneous  surface  rainfall  is  represented  by  the  contour  lines.  The 
plotted  domain  is  Qpi0t  =  180  x  120  x  17.5  km3.  (In  Matlab:  azimuth  =  -110,  elevation  =  10). 
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Figure  8:  FM06:  perspective  view  of  the  squall  line  (grey  scale)  at  1800,  3600,  and  7200  seconds.  Stabilization 
achieved  by  VMS+DC  (left  column)  and  by  HV  with  v 4  =  1  x  108m4s-1  (right  column).  The  velocity  vectors  at 
the  surface  are  plotted  at  every  7  grid  points.  The  instantaneous  surface  rainfall  is  represented  by  the  contour  lines. 
The  plotted  domain  is  Qpiot.  =  300  x  60  x  17.5  km3.  (In  Matlab:  azimuth  =  -110,  elevation  =  10). 


17 


the  same  rate  of  decay  and  extremely  close  values  along  the  whole  18000  s.  The  absolute  values 
and  the  trends  are  congruent  with  the  literature2. 

Looking  at  the  curves  of  min(gc)  and  min(gr)  in  Fig.  10,  we  see  that  our  method  is  consistently 
controlling  the  undershoots  of  the  tracers.  It  is  interesting,  however,  to  notice  how  the  gap  between 
VMS+DC  and  any  HV  for  the  current  test  is  no  longer  as  pronounced  as  it  was  for  Squall  line  I 
(Fig.  4).  We  may  have  found  the  proper  tuning  of  the  HV  coefficients  that  greatly  limit  min(gr), 
although  it  still  does  not  limit  min(gc). 

The  line-averaged  contours  of  the  cloud  content  are  plotted  in  Fig.  11.  As  expected  by  looking 
at  the  time  series,  the  VMS+DC  and  HV  solutions  are  similar,  with  a  relatively  small  difference 
in  the  value  of  qc.  The  overall  shape  (characteristic  anvil),  the  height  of  the  lifting  condensation 
level,  and  the  vertical  extension  of  qc  agree  with  Weisman  et  al.  (1988)  and  Weisman  and  Rotunno 
(2004)  for  both  solutions.  The  amount  of  accumulated  rain  confirms  the  similarities.  We  see  this 
in  Fig.  12  where  the  contours  of  total  accumulated  rainfall  are  plotted  at  tf  =  18000  s.  We  notice 
that  the  structure  of  the  precipitation  resolved  by  our  method  does  not  resemble  the  shape  of  the 
squall  line  as  much  as  it  did  in  Fig.  5.  Nevertheless,  the  irregular  pattern  is  supported  by  Frame 
and  Markowski  (2006),  although  with  a  difference  of  a  few  millimeters  in  the  maximum  value  of 
precipitated  water  in  both  VMS+DC  and  HVa. 

4.  Conclusions 

We  have  presented  a  dynamic  and  parameter-free  alternative  to  the  classical  fourth  order  hyper¬ 
viscosity  to  stabilize  the  solution  of  the  coupled  transport  equations  of  water  tracers  to  simulate  3D 
deep  convection  using  spectral  elements.  The  current  method,  based  on  the  Variational  Multiscale 
Stabilization  (VMS)  with  crosswind  discontinuity  capturing  (DC),  has  shown  to  have  certain  im¬ 
portant  properties  that  make  it  a  suitable  stabilizing  alternative  to  hyper- viscosity  in  the  simulation 
of  atmospheric  problems  modeled  by  a  system  of  coupled  transport  equations  with  phase  change. 
Most  importantly,  VMS+DC  is  parameter-free  -a  clear  advantage  from  the  point  of  view  of  the 
user  of  the  code;  it  is  numerically  consistent;  and  it  is  anisotropic.  Furthermore,  it  improves  the 
positivity  preserving  properties  of  the  spectral  element  solution  of  the  advection  problem  without 
the  need  for  additional  filters.  In  the  simulation  of  moist  convection,  this  property  is  of  major 
importance  although  it  remains  a  challenge  to  construct  positivity  preserving  schemes  at  higher 
order  (e.g.,  see  Marras  et  al.  (2012)  for  one  possibility  of  tackling  this  problem).  Although  our 
approach  is  not  fully  monotonicity  preserving  (at  least  not  for  spectral  element  basis  functions 
beyond  second  order),  it  still  helps  in  minimizing  the  undershoots.  Obtaining  negative  tracers  is 
of  major  concern  in  weather  simulations  and  climate  modeling.  With  respect  to  parallel  efficiency, 
because  each  Laplacian  requires  one  communication,  fewer  communications  are  required  by  VMS 


2  Because  Frame  and  Markowski  (2006)  do  not  provide  this  information,  we  rely  on  a  comparison  with  Weisman 
et  al.  (1988)  and  Weisman  and  Rotunno  (2004)  whose  initial  state  is  the  same. 
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Figure  9:  FM06:  Time  evolution  of  max(w)  (top)  and  min(w)  (bottom). 
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Figure  10:  FM06:  time  evolution  of  min(qc)  (top)  and  min(gr)  (bottom) 
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Figure  11:  FM06:  line  averaged  vertical  cross  sections  of  qc  at  1800,  3600,  and  7200  seconds.  Stabilization  achieved 
by  VMS+DC  (left  column)  and  by  HV  with  i/4  =  1  x  10s  m4s^1  (right  column). 
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Figure  12:  FM06:  total  accumulated  rainfall  (mm)  at  tf  =  18000s.  Left:  VMS+DC.  Right:  HVa  with  ux  =  vy  = 
1  x  109  m4  s_1,  vz  =  3  x  105  m4  s_1. 


Table  2:  Salient  properties  of  VMS+DC  and  HV 


Method 

Parameter- free 

Anisotropic 

Consistent 

Monotonicity 

Parallel  communication 

VMS  +  DC 

V 

V 

V 

Improved 

As  V2 

HV 

X 

X 

As  V2q 

and  DC  than  by  high  order  dissipative  operators.  As  we  try  to  approach  the  exascale  range  of 
computing,  fewer  communications  are  essential  for  faster  and  more  scalable  codes.  We  summarize 
the  most  salient  properties  of  VMS+DC  in  Table  2. 

Although  the  method  hereby  described  was  originally  designed  around  a  second  order  operator, 
we  do  not  foresee  a  reason  why,  by  a  dimensionally  consistent  correction,  it  could  not  be  adapted  for 
higher  order  operators  in  case  one  were  interested  in  extending  this  machinery  to  a  parameter-free 
higher  order  hyper-viscosity-like  dissipation.  In  this  case,  however,  the  cost  of  parallel  communica¬ 
tion  would  increase. 

Finally,  although  it  was  built  for  finite  and  spectral  elements,  the  discontinuity  capturing  dis¬ 
sipation  DC  could  be  constructed,  with  little  effort,  for  alternative  numerical  techniques  such  as 
finite  volume  and  finite  difference  methods. 
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